##########################################################################################

library(dplyr)
library(ggplot2)
library(data.table)
library(RColorBrewer)
library(optparse)
library(ggpubr)

##########################################################################################

option_list <- list(
    make_option(c("--gene"), type = "character") ,
    make_option(c("--mut_rate_gene_file"), type = "character") ,
    make_option(c("--mut_rate_point_file"), type = "character") ,
    make_option(c("--images_path"), type = "character")
)

if(1!=1){
    
    gene <- "TP53"
    work_dir <- "~/20220915_gastric_multiple/dna_combine_20221213/"
    mut_rate_gene_file <- paste(work_dir,"/images/mutRate/MutRate.tsv",sep="")
    mut_rate_point_file <- paste(work_dir,"/images/mutRate/MutRate.RecurrentPoint.tsv",sep="")
    images_path <- paste(work_dir,"/images/mutRatePlot",sep="")

}

###########################################################################################

parseobj <- OptionParser(option_list=option_list, usage = "usage: Rscript %prog [options]")
opt <- parse_args(parseobj)
print(opt)

gene <- opt$gene
mut_rate_gene_file <- opt$mut_rate_gene_file
mut_rate_point_file <- opt$mut_rate_point_file
images_path <- opt$images_path

dir.create(images_path , recursive = T)

###########################################################################################

col <- c(
  brewer.pal(9,"YlGnBu")[6],
  rgb(234,106,79,alpha=255,maxColorValue=255),
  rgb(203,24,30,alpha=255,maxColorValue=255),
  rgb(255,0,0,alpha=255,maxColorValue=255)
  )

names(col) <- c("IM" , "IGC" , "DGC" , "GC")

###########################################################################################

dat_mutRateGene <- data.frame(fread( mut_rate_gene_file ))
dat_mutRatePoint <- data.frame(fread( mut_rate_point_file ))

###########################################################################################

dat_mutRateGene <- subset(dat_mutRateGene , Hugo_Symbol==gene)
dat_mutRatePoint <- subset(dat_mutRatePoint , Hugo_Symbol==gene)
dat_mutRatePoint$Hugo_Symbol <- dat_mutRatePoint$vid
dat_mutRatePoint <- dat_mutRatePoint[,-3]

dat_plot <- rbind( dat_mutRateGene , dat_mutRatePoint )
dat_plot$Class <- factor( dat_plot$Class , levels = c("IM" , "GC" , "IGC" , "DGC") , order = T )
dat_plot$value_text <- paste0( round(dat_plot$MutRate , 2) * 100 , "%") 

###########################################################################################

trans <- function(num){
    up <- floor(log10(num))
    down <- round(num*10^(-up),2)
    text <- paste("p == ",down," %*% 10","^",up)
    return(text)
}

###########################################################################################
## 计算P值
dat_plot$p.value = ""
dat_plot$p_text = ""
dat_plot <- subset( dat_plot , Hugo_Symbol %in% unique(dat_plot[dat_plot$MutNum > 2,"Hugo_Symbol"]))

class_compare <- c("IM" , "IGC" , "DGC")

result <- c()
for(geneN in unique(dat_plot$Hugo_Symbol)){

    print(geneN)

    for( i in 1:(length(class_compare)-1) ){
        class1 <- class_compare[i]
        for( j in (i+1):length(class_compare) ){
            class2 <- class_compare[j]
            tmp_1 <- subset( dat_plot , Hugo_Symbol == geneN & Class %in% class1 )
            tmp_2 <- subset( dat_plot , Hugo_Symbol == geneN & Class %in% class2 )

            if(nrow(tmp_1)==0){
                tmp_1 <- tmp_2
                tmp_1$MutNum <- 0
                tmp_1$MutRate <- 0
                tmp_1$value_text <- ""
            }

            if(nrow(tmp_2)==0){
                tmp_2 <- tmp_1
                tmp_2$MutNum <- 0
                tmp_2$MutRate <- 0
                tmp_2$value_text <- ""
            }

            tmp <- rbind( tmp_1 , tmp_2 )
            tmp_fisher <- matrix(c(tmp$MutNum , tmp$SampleNum - tmp$MutNum) , ncol = 2)
            p <- fisher.test(tmp_fisher)$p.value

            if( p < 0.001 ){
                p_text <- trans(p)
            }else{
                p_text <- paste0( "p == " , round(as.numeric(p) , 3) ) 
            }

            tmp$p.value <- p
            tmp$p_text <- p_text

            result <- rbind( result , tmp )
        }
    }
}

result$Hugo_Symbol <- factor( result$Hugo_Symbol , 
    levels = unique(result[order(result$MutRate , decreasing=T),"Hugo_Symbol"]) , order = T)
result$Class <- factor( result$Class , levels = c("IM" , "IGC" , "DGC") , order = T )
out_name <- paste0( images_path , "/MutRate_" , gene , ".IM_IGC_DGC.tsv" )
write.table( result , out_name , row.names = F , sep = "\t" , quote = F )

###########################################################################################

dat_use <- unique( result[,1:6] )

plot <- ggplot( data = dat_use , aes( x = Class , y = MutRate , fill = Class ))+
geom_bar(position = "stack", stat = "identity") + 
theme_bw()+
labs(x="",y="Mutation Rate")+
facet_grid(.~Hugo_Symbol) +
theme(panel.grid = element_blank())+
scale_fill_manual(values=col) +
ylim(0,1)+
geom_text(aes(label=value_text) , position=position_stack(vjust = 0.5) , size=2.5 , color="black")+
theme(panel.background = element_blank(),#设置背影为白色#清除网格线
            legend.position ='right',
            legend.title = element_blank() ,
            panel.grid.major=element_line(colour=NA),
            plot.title = element_text(size = 12,color="black",face='bold'),
            legend.text = element_text(size = 12,color="black",face='bold'),
            axis.text.y = element_text(size = 12,color="black",face='bold'),
            axis.title.x = element_text(size = 12,color="black",face='bold'),
            axis.title.y = element_text(size = 12,color="black",face='bold'),
            strip.text.x = element_text(size = 7 , face = 'bold'),
            axis.ticks.x = element_blank(),
            axis.text.x = element_text(size = 6,color="black",face='bold') ,
            axis.line = element_line(size = 0.5))

out_name <- paste0( images_path , "/MutRate_" , gene , ".IM_IGC_DGC.pdf" )
ggsave( out_name , plot , width = (3 + 1.5 * length(unique(dat_plot$Hugo_Symbol))) , height = 5 )

###########################################################################################
## 画标记p值的图
barplot_pval <- function(genename,data,out_name){
  ## 这里提取了一下需要的内容 至于选择的顺序不知道这样的处理是否合理

  testdata <- data[which(data$Hugo_Symbol==genename),]
  barplotdata <- testdata[order(testdata$Class,decreasing = T),]
  barplotdata <- barplotdata[c(1,3,5),c(1,5,6)]
  pvaldata <- testdata[c(1,3,5),c(1,7)]
  
  ##设计横杠的位置
  if (max(barplotdata$MutRate)>0.9){
    maxheight <- 1.1
    UNIT <- 0.1
  }else {
    maxheight <- round(max(barplotdata$MutRate),3)
    spare <- (1-maxheight)
    UNIT <- round(spare/3,3)
  }

  height1 <- maxheight + UNIT
  height2 <- height1 + UNIT
  height3 <- height2 + UNIT
  line <- data.frame(x=c(0.5,0.5,3.5),xend=c(2,3.5,2),
                     y=c(height1,height3,height2),yend=c(height1,height3,height2))
  
  ##p值  字的位置和颜色
  trans <- function(num){
    if(num < 0.001){
        up <- floor(log10(num))
        down <- round(num*10^(-up),2)
        text <- paste("p == ",down," %*% 10","^",up)
    }else{
        text <- paste("p == ",round(num , 3))
    }
    return(text)
  } #最后的效果是在类似expression()中呈现的 "=="表示等于 "%*%"表示✖
  
  TEXTUNIT <- UNIT/1.7
  
  if((pvaldata$p.value[1])<0.05){
    curcol="#CA171E"
  }else{
    curcol="#000000"
  }
  pval12 <- data.frame(x=1.2,y=height1+TEXTUNIT,label=trans(pvaldata$p.value[1]),col=curcol)
  
  if((pvaldata$p.value[2])<0.05){
    curcol="#CA171E"
  }else{
    curcol="#000000"
  }
  pval13 <- data.frame(x=2,y=height3+TEXTUNIT,label=trans(pvaldata$p.value[2]),col=curcol)
  
  if((pvaldata$p.value[3])<0.05){
    curcol="#CA171E"
  }else{
    curcol="#000000"
  }
  pval23 <- data.frame(x=2.8,y=height2+TEXTUNIT,label=trans(pvaldata$p.value[3]),col=curcol)
  
  ##设计竖杠的位置
  UNIT <- UNIT/2
  group12 <- data.frame(x=0.5,xend=0.5,x2=2,xend2=2,y=height1,yend=(height1-UNIT))
  group13 <- data.frame(x=0.5,xend=0.5,x2=3.5,xend2=3.5,y=height3,yend=(height3-UNIT))
  group23 <- data.frame(x=2,xend=2,x2=3.5,xend2=3.5,y=height2,yend=(height2-UNIT))
  
  plot <- ggbarplot(barplotdata, x="Class",y="MutRate",
            label = barplotdata$value_text, fill = "Class" ,
            lab.pos="in") +
    xlab(label='') +
    scale_fill_manual(values=col) +
    ylab(label="Mutation Rate") +
    ggtitle(genename)  +
    geom_segment(aes(x = x, y = y, xend = xend, yend = yend),data=line) +
    geom_segment(aes(x = x, y = y, xend = xend, yend = yend),data=group12) +
    geom_segment(aes(x = x2, y = y, xend = xend2, yend = yend),data=group12) + 
    geom_segment(aes(x = x, y = y, xend = xend, yend = yend),data=group13) +
    geom_segment(aes(x = x2, y = y, xend = xend2, yend = yend),data=group13) +
    geom_segment(aes(x = x, y = y, xend = xend, yend = yend),data=group23) +
    geom_segment(aes(x = x2, y = y, xend = xend2, yend = yend),data=group23) +
    geom_text(aes(x=x,y=y,label=label),data=pval12,size=4.2,parse = TRUE,colour = pval12$col) +
    geom_text(aes(x=x,y=y,label=label),data=pval13,size=4.2,parse = TRUE,colour = pval13$col) +
    geom_text(aes(x=x,y=y,label=label),data=pval23,size=4.2,parse = TRUE,colour = pval23$col) +
    theme(panel.background = element_blank(),#设置背影为白色#清除网格线
            legend.position ='none',
            legend.title = element_blank() ,
            panel.grid.major=element_line(colour=NA),
            plot.title = element_text(size = 12,color="black",face='bold'),
            legend.text = element_text(size = 12,color="black",face='bold'),
            axis.text.y = element_text(size = 7, color="black",face='bold'),
            axis.title.x = element_text(size = 12,color="black",face='bold'),
            axis.title.y = element_text(size = 12,color="black",face='bold'),
            strip.text.x = element_text(size = 7 , face = 'bold'),
            axis.ticks.x = element_blank(),
            axis.text.x = element_text(size = 10,color="black",face='bold') ,
            axis.line = element_line(size = 0.5))

     
  ggsave( out_name , plot , width = 4 , height = 5 )
}


metadata <- result #load data

for( genename in unique(metadata$Hugo_Symbol) ){

    out_name <- paste0( images_path , "/MutRate_" , genename , ".IM_IGC_DGC.addPvalue.pdf" )
    barplot_pval(genename,metadata,out_name)

}

